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Calculating Tumor Trajectory and Dose-of-the-day Using Cone-Beam CT 
Projections 
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Department of Radiation Oncology, University of Colorado School of Medicine 

Purpose: Cone-beam CT (CBCT) projection images provide anatomical data in real-time over several 
respiratory cycles, forming a comprehensive picture of tumor movement. We developed and validated a 
method which uses these projections to determine the trajectory of and dose to highly mobile tumors during 
each fraction of treatment. 

Methods: CBCT images of a respiration phantom were acquired, the trajectory of which mimicked a 
lung tumor with high amplitude (up to 2.5 cm) and hysteresis. A template-matching algorithm was used 
to identify the location of a steel BB in each CBCT projection, and a Gaussian probability density function 
for the absolute BB position was calculated which best ht the observed trajectory of the BB in the imager 
geometry. Two modihcations of the trajectory reconstruction were investigated: first, using respiratory phase 
information to rehne the trajectory estimation (Phase), and second, using the Monte Carlo (MC) method to 
sample the estimated Gaussian tumor position distribution. The accuracies of the proposed methods were 
evaluated by comparing the known and calculated BB trajectories in phantom-simulated clinical scenarios 
using abdominal tumor volumes. 

Results: With all methods, the mean position of the BB was determined with accuracy better than 0.1 
mm, and root-mean-square (RMS) trajectory errors averaged 3.8±1.1% of the marker amplitude. Dosimetric 
calculations using Phase methods were more accurate, with mean absolute error less than 0.5%, and with 
error less than 1% in the highest-noise trajectory. MC-based trajectories prevent the over-estimation of dose, 
but when viewed in an absolute sense, add a small amount of dosimetric error (<0.1%). 

Conclusions: Marker trajectory and target dose-of-the-day were accurately calculated using CBCT pro¬ 
jections. This technique provides a method to evaluate highly-mobile tumors using ordinary CBCT data, and 
could facilitate better strategies to mitigate or compensate for motion during SBRT. 

This manuscript was submitted to Medical Physics 


I. INTRODUCTION 

As radiotherapy becomes more conformal, the nega¬ 
tive effects of tumor motion become more pronounced^. 
Highly targeted techniques, such as Stereotactic Body 
Radiotherapy (SBRT), rely on the establishment of a 
very steep dose gradient around the tumor in order to 
reduce normal tissue effects^. Motion blurs the bound¬ 
ary between tumor and normal tissue, making this dose 
gradient harder to establish. Interfraction motion can 
largely be addressed by image-guided therapy techniques 
such as cone-beam CT (CBCT)^; however, intrafraction 
motion presents a more complicated problem. Specifi¬ 
cally, motion makes it more difficult to ensure that the 
target receives the full prescribed dose, and also adds un¬ 
certainty when calculating the dose-of-the-day based on 
pre-treatment imaging. 

A common approach to account for regular intrafrac¬ 
tion motion is to use 4-Dimensional CT (4DCT) imaging 
to characterize the full range of tumor motion, and to 
treat a target volume which encompasses this range^’'*. 
This guarantees full coverage of the tumor but increases 
dose to normal tissue in the tumor vicinity. Another ap¬ 
proach is to monitor the tumor position (either directly or 
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using surrogates), and to treat only when the tumor lies 
within some specified region (i.e. gating)®’’’. More recent 
advances involve modifying or re-directing the treatment 
beam in real-time in order to compensate for intrafrac¬ 
tion motion®’®. 

More accurate methods to account for tumor motion 
are based on knowing the instantaneous tumor location 
over time (i.e. the tumor trajectory). Several meth¬ 
ods exist to determine this. If the tumor is implanted 
with radiopaque hducial markers, the position can be tri¬ 
angulated using stereoscopic imaging from two or more 
view angles’®*’®. Similarly, the tumor can be implanted 
with electromagnetic transponders’’’”, or can be moni¬ 
tored in the treatment suite using MRI images’®. One 
major drawback of these methods is that they rely on 
technologies or treatments not commonly available to the 
majority of radiotherapy centers, or are only available for 
use in certain tumor sites. 

CBCT is a widely used pre-treatment imaging modal¬ 
ity which captures a relatively low-quality 3D image of 
the patient anatomy in the treatment position®. Typi¬ 
cal clinical CBCT systems capture roughly 650 projec¬ 
tion images, which are used to reconstruct the 3D image. 
These projections are typically acquired over a 30-120 
second timescale, and the reconstructed images exhibit 
significant blurring due to respiratory motion. However, 
the projection images themselves provide a comprehen¬ 
sive picture of the anatomical movement, as they typ¬ 
ically observe the motion of several respiratory cycles. 
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For tumors implanted with fiducial markers, these pro¬ 
jections can reveal the motion of the target volume dur¬ 
ing treatment. Unfortunately, the clinical utility of these 
projections is often limited as these data are incomplete 
with respect to the tumor position, since the process of 
projection condenses the location of an object in 3D space 
(x,y,z) to only two dimensions (u,v) in the imager ge¬ 
ometry. Calculation of the tumor trajectory from these 
projections would allow for determination of tumor mo¬ 
tion and dose-of-the-day estimation for tumors implanted 
with fiducial markers, which is often the case for tumors 
of the lung, liver, or pancreas^®"^®. Additionally, this 
information could be used in the context of Adaptive 
Radiotherapy”^^'®^. 


The goal of this work was to calculate the tumor tra¬ 
jectory using CBCT projections, and to validate the ac¬ 
curacy of this trajectory in calculating the tumor dose. 
Mathematically, determining the instantaneous 3D po¬ 
sition of the tumor using the 2D projections acquired 
during CBCT is an under-defined problem; in each pro¬ 
jection we have two measured data points {u,v position 
within the imager geometry) yet wish to calculate three 
quantities {x,y,z position of the tumor). The projection 
measurement confines the position to lie along a line be¬ 
tween the imaging source and the point of detection, but 
cannot by itself determine the depth of the position along 
this line. Methods have been developed to estimate the 
tumor position using orbiting CBCT projections, mostly 
in the context of tracking/targeting mobile tumors in 
real-time^^'^^. These methods choose the tumor location 
by maximizing the expectation value of tumor position 
(which maximizes the probability of hitting the tumor). 
However, this may not be appropriate for dosimetric eval¬ 
uation, since the most likely tumor location is generally 
in the center of the planning target volume (PTV) and 
may lead to an over-estimation of the tumor dose. More¬ 
over, tumors which undergo large respiratory motion of¬ 
ten exhibit a cyclical trajectory. Previous methods are 
built around the assumption that tumor position is dis¬ 
tributed randomly about the daily mean location, leading 
to a Gaussian distribution of random position errors^®. 
However, this assumption is broken in the presence of 
non-random respiratory motion, and a Gaussian distri¬ 
bution of these positions will contain relatively large un¬ 
certainty. For cyclical motion, the true uncertainty is 
much smaller, since there is a large correlation between 
position and respiratory phase. 


In this work, we extended these methods for more ac¬ 
curate calculation of tumor dosimetry in highly-mobile 
tumors by incorporating information regarding the res¬ 
piratory phase. Using a motion phantom, we validated 
the accuracy of these methods, and quantified the uncer¬ 
tainty in using GBCT projections to calculate dose-of- 
the-day in mobile tumors. 


II. METHODS 

To quantify the motion of the target, the location of a 
marker in 3D phantom geometry is calculated by exam¬ 
ining the position of that marker in a series of orbiting 
2D projection images. Our method builds on the work of 
Poulsen et and proceeds as follows. First, we acquire 
CBCT projections of a phantom in which the position of 
a small fiducial marker is visible. Next, the position of 
the fiducial marker in these projections is determined au¬ 
tomatically using a template-matching algorithm. Then, 
the acquired images are binned into different phases of 
the respiratory cycle. An estimate of the likely distri¬ 
bution of marker positions in each phase is built using 
a Gaussian approximation of uncertainty^^’®®. Finally, 
this Gaussian distribution is used to estimate the unre¬ 
solved component of position, either by calculating the 
most likely position or by using Monte Carlo methods to 
sample this distribution. 


A. Definition of geometry 

Let (x,y,z) represent an orthogonal basis in 3D space, 
with the origin corresponding to the treatment and imag¬ 
ing isocenter. In this work, we consider the patient to 
lie in the head-first supine position, and choose (x,y,z) to 
correspond to the anterior-posterior (AP) or vertical, left- 
right (LR) or lateral, and superior-inferior (SI) or longi¬ 
tudinal directions, respectively. The kV source/imager is 
oriented orthogonal to the z-axis, and orbits the patient 
in the x/y plane (rotation about the z-axis). 9 represents 
the angle of the kV imaging source with respect to the x 
axis. Let (u,v) describe the axes of the 2D projection im¬ 
age, where u is parallel to the x/y plane and v is parallel 
to the z-axis. Note that this coordinate system matches 
those typically described in the radiotherapy literature 
(e.g. the seminal work of Feldkamp et aF/] however, 
particular choice of the relative directions of x, y, z, u, v, 
and 9 may alter the relationship and/or sign dependence 
of these quantities. 


B. Projection operator 

Let P define the projection of a position (x,y,z) into the 
(u,v) plane at an angle 9; in other words, P represents 
the acquisition of a CBCT projection. SDD is the source- 
to-detector distance, p is the image pixel size, SAD is 
the source-to-axis distance, and o„,Ot, are the projected 
locations of the rotational isocenter in the u,v coordinates 
(the so-called piercing-point or principal point). In Eqs 
1-3, u and v are given in units of pixels in imaging space. 


{u,v) = P{x,y,z,9) (1) 
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SDD X s\n{0) + y cos{0) 
p S'^-D — [a; cos(0) — ysin(0)] 


K = 



-(p+e/i-ro)^B(p+e;x-ro)/2 


( 7 ) 


p — [a;cos(0) — j/sin(0)] “ 

Consider the following clinical scenario. A patient is 
implanted with several radiopaque fiducial markers in or¬ 
der to assist in the delineation and localization of this 
tumor for treatment. A CBCT is acquired in order to vi¬ 
sualize the patient anatomy on the day of treatment, and 
involves the acquisition of several hundred images which 
contain the projection of that anatomy into 2D. The ac¬ 
quisition of this CBCT scan is represented by performing 
P(0) for values of 0 from 0° to 360°, and for each fidu¬ 
cial marker, each image contains a measurement of the 
marker position ug and vg. 

C. Marker position 


g{t) = Ke-(‘-™)V( 2 <T)^ (8) 

Here, e describes a unit vector pointing from the ob¬ 
served marker position to the imaging source, and p de¬ 
scribes the point on the imager at which the marker was 
detected (converted to patient geometry). 

Using the observed marker positions {ui, Vi) in each im¬ 
age i, optimization was performed to estimate the mean 
position and covariance matrix which describes the prob¬ 
abilistic marker positions. By integrating Eq. 8, it can 
be shown that the probability of finding the marker at 
a given detector position is linearly proportional to the 
product K ■ a. The optimization finds a mean position 
and covariance matrix which maximizes the total prob¬ 
ability of finding the markers in positions observed in 
each image i. The optimization was performed using the 
Nelder-Mead simplex search method^’’. 


We used the formalism of Poulsen et to construct a 
3D Gaussian probability density function (PDF) describ¬ 
ing the marker position. A full derivation of this Gaus¬ 
sian formalism is given in Ref 22, and for convenience of 
the reader, we use the same notation when possible. This 
method uses maximum-likelihood optimization to calcu¬ 
late a PDF which best describes the observed motion of 
a marker in projection images. Then, using this PDF, 
the unresolved component of marker position is deter¬ 
mined by computing the expectation value of this PDF 
along the line between the x-ray source and the detection 
point within the image (Fig 1). 

A 3D Gaussian PDF is defined by 9 quantities: mean 
location (rg) = (xg, yo)-^o): variance (a^, cry, and co- 
variance (a^y, CTyz, (Jzx)- The variance and covariance to¬ 
gether form the covariance matrix A. From this, the 3D 
PDF for marker position (r) = {xo,yo,zo) was given by: 


G{x,y,z) 


det B rp_x-)TB(r-r„)/2 

( 271)3 


( 4 ) 


Here, B denotes the inverse covariance matrix B = 
=-i 

A , and T denotes the transpose operator. The func¬ 
tional form of Eq. 4 along the line from the imaging 
source to detection point (defined here as g) was de¬ 
scribed by the following: 


cr = 


e^Be 


( 5 ) 


(n) - pj'^Be 

eABe 


( 6 ) 


[rg, B] = arginin^ - \og(Kiai) (9) 

ro,S i 


D. Respiratory phase 

Respiratory phase was determined using the SI dis¬ 
placement of the fiducial markers. The SI position is 
well-determined in the GBCT projections for locations 
near the rotational isocenter, since the only error derives 
from magnification effects related to (x, y) errors (Eq. 3), 
which are typically on the order of 1% or less. This in¬ 
ternal surrogate information precluded the need to define 
respiratory phase based on some other surrogate, such as 
chest markers or diaphragm movement. 

First, the marker trajectory in the v direction was 
smoothed using a five-sample moving average, and dis¬ 
crete velocity was calculated by taking the difference be¬ 
tween smoothed positions in sequential projections. It 
is assumed that these marker positions represent tumor 
trajectory in a head-first supine orientation, with the -l-v 
axis pointing in the superior direction. End-exhale (fee) 
was identified from projections where the marker velocity 
switched from positive to negative; likewise, end-inhale 
(vei) was identified where velocity went from negative to 
positive. 

Projections were binned into one of 10 respiratory 
phases based on the velocity and location of the marker 
relative to the end-exhale and end-inhale positions. 
Phases 1-5 (inhale) consisted of projections with nega¬ 
tive velocity, grouped into uniformly spaced bins between 
subsequent Vee and Vei coordinates. Phases 6-10 (exhale) 
consisted of projections with positive velocity, grouped 
into linearly spaced bins between subsequent Vei and Vee 
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FIG. 1. Determination of tumor position using 3D Gaussian probability. A) Using all 2D projections, a 3D Gaussian distribution 
is found which best fits the observed marker trajectory. For each measured position (i.e. each projection image), the position is 
determined by finding the expectation of this Gaussian distribution along the line between the x-ray source and the detection 
point. B) In the phase-resolved method, a 3D distribution is independently calculated for each phase of the respiratory cycle. 


coordinates. Projections at the beginning (or end) of the 
motion trace (i.e. those that did not contain measure¬ 
ments of the full cycle) used Vei (or Vee) coordinates of 
the nearest full phase. 


E. Position determination 


Four methods were used to determine the 3D marker 
trajectory from the projection images. These meth¬ 
ods were implemented using in-house software developed 
with Matlab (The Mathworks, Natick MA). In Method 
1 (Base), position estimation followed the method of 
Poulsen et aP^. Data from all projections was used dur¬ 
ing optimization to calculate a 3D Gaussian distribution 
which described the probability of finding the marker in 
a given 3D location (Eq. 9). This equation was parame¬ 
terized along the detection ray (the line between the ra¬ 
diation source and the imager detection point), and the 
marker position along this line was given hy t = fx, where 
fx equals the expectation value of the distance along this 
line (Eqs. 5-8). 

In Method 2 (MC), Monte Carlo methods were used 
to select the marker position in 3D space in place of the 
expectation value /i. First, the detection ray was dis¬ 
cretized into 10,000 points, spaced linearly between the 
source and detector. The PDF of the ID Gaussian (Eq. 
8) was integrated across the discretized ray to form a 
normalized cumulative probability distribution function. 


CDF{t^) = 




( 10 ) 


EZXgiu) 

Next, a random number ^ was uniformly chosen be¬ 


tween 0 and 1. The position of the marker along the ray 
was chosen to be the point T which satisfied CDF{T) = 
and this solution was found numerically. 

To increase the accuracy of tracking for targets with 
high respiratory motion, the trajectory of the marker 
was calculated using a phase-resolved Gaussian approxi¬ 
mation (Method 3, Phase). All projections were binned 
into one of 10 respiratory phases, and Eq. 9 (optimiza¬ 
tion of Gaussian distribution) was applied to data from 
each phase independently. This resulted in a set of 10 
phase-specific mean positions and inverse covariance 

matrices Bp, or in other words, 10 distinct distributions 
describing the position of the marker during each phase 
of respiration (Fig IB). In Method 4 (PhaseMG), both 
the Phase and MG methods were applied simultaneously. 


F. Motion phantom 

In order to validate the proposed methods, a phantom 
was used to simulate respiratory tumor motion (Quasar 
Phantom, Modus Medical Devices, Inc.). This phantom 
contains a 1-mm-diameter steel BB, the trajectory of 
which is shown in Fig 2 for various superior-inferior (Sup- 
Inf) amplitudes of motion. The trajectory is similar to 
that of a lung tumor trajectory with high amplitude and 
hysteresis". A GBCT scan of this phantom was acquired 
using the imaging capabilities of a Varian Truebeam ac¬ 
celerator (Varian Medical Systems, Palo Alto, GA). 656 
projection images were acquired using a half-fan trajec¬ 
tory {9 = —180° —)■ 180°). The phantom was set to a 
rate of 11 breaths-per-minute, with a maximum range 
in the SI direction of 2.5 cm (resulting in maximum AP 
and LR ranges of 1.2 and 0.6 cm, respectively). Multiple 
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FIG. 2. Trajectory of Motion Phantom. The motion of the 1-mm-diameter steel BB embedded within the motion phantom 
is shown relative to the isocenter of the CBCT imaging system. Trajectories are shown for amplitudes of motion in the Sup- 
Inf direction ranging from 5 25 mm. For clarity, the trajectories are also shown projected onto the XY, XZ, and YZ axes. 
Trajectories are shown relative to a standard head-first supine orientation (i.e. the primary axis of motion in the Sup-Inf 
direction). 


projection image datasets were acquired with trajectory 
amplitudes ranging in 1 mm increments from 1 mm to 
2.5 cm. 

In each image, the location of the BB was determined 
automatically using a template matching algorithm^®. A 
discretized 2D Gaussian kernel was used as a template 
image, with width equal to the radius of the BB. In each 
image, the normalized cross correlation was computed 
between the projection image and the template^®. The 
location of the BB in the image was given by the location 
where the normalized cross correlation was maximized. 
In cases where the template matching algorithm failed, 
the position was determined manually. 

Using the positions of the BB in the projection images, 
the trajectory of the BB in 3D space was calculated using 
methods 1-4 (Sect 2.5). The accuracy of these methods 
was calculated by comparing the calculated trajectory to 
the true trajectory (shown in Fig 2). To simulate more 
difficult and unpredictable clinical scenarios, additional 
(simulated) trajectories were created. Gaussian noise was 
added to the true marker trajectory by adding normally 
distributed random values to the x, y, and z locations at 
each point in time. Simulated trajectories were created 
by varying the standard deviation of this noise from 0 
to 1 mm in 0.1 mm increments between datasets. These 
altered datasets were projected into 2D image space (Eqs. 
1-3), and Methods 1-4 were used to calculate the marker 
trajectory as described previously. 


G. Dosimetric model 

The goal of this study was to evaluate the accuracy 
of dosimetric calculations based on reconstructed trajec¬ 
tories. To that end, the dosimetry of the reconstructed 
trajectories was compared to the dosimetry of the true 
trajectory. To evaluate the dosimetry of the true and 
reconstructed phantom motion trajectories, a phantom- 
simulated clinical scenario was considered. A clinical tar¬ 
get volume (CTV) was assumed to undergo the motion 
trajectories under analysis. Treatment was prescribed 
to a planning target volume (PTV), which included the 
CTV plus an isotropic margin. To simplify the model, 
and to exaggerate differences between different methods, 
the dose within the PTV was assumed to equal the pre¬ 
scription, with zero dose outside the PTV. The PTV (and 
resulting dose distribution) was stationary over time, 
while the CTV was translated within the PTV based on 
either the true phantom trajectory or the reconstructed 
trajectories. For each of the 656 points in the trajectory 
(corresponding to the 656 projection images), the posi¬ 
tion of the CTV relative to its center was given by the 
instantaneous trajectory position relative to the mean 
position. Based on this position, the dose was calculated 
for each voxel in the CTV (i.e. based on the binary dose 
model, or in other words assigning l/656th of the pre¬ 
scription dose for voxels inside the PTV, and zero for 
those outside). Independent simulations were performed 
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Trajectory Amplitude (mm) Trajectory Amplitude (mm) Trajectory Error (cm) 


FIG. 3. Errors in calculating trajectory using orbiting views. Results were calculated with no corrections (Base), with phase 
corrections (Phase), with Monte Carlo sampling (MC) and with both phase and Monte Carlo (PhaseMC). Results are shown 
for phantom motion amplitudes in the SI direction ranging from 1 mm to 25 mm. A) Root-mean-squared (RMS) errors in 
the trajectory reconstruction using different methods. Phase and PhaseMC methods reduce RMS errors by half. B) Relative 
fraction of trajectory errors greater than 3 mm. Phase methods reduce the fraction of large errors, and Monte Carlo methods 
increase the percentage of large errors slightly. C) Histogram of trajectory errors in all projections across all acquired datasets. 
Phase methods demonstrate increased accuracy. 


for CTV-to-PTV margins ranging from 0 mm to 15 mm, 
and were deliberately set less than the amplitude of mo¬ 
tion so as to allow for differences between the proposed 
methods to become apparent. The goal of our analysis 
was to quantify the accuracy of dosimetric calculations 
using CBCT-based trajectory reconstruction. To that 
end, our binary dose model serves as a conservative esti¬ 
mate of the inaccuracies of this approach, as any errors in 
position determination will be reflected in the idealized 
nature of the target dosimetry simulation. 


H. Validation metrics 

We analyzed the accuracy of the proposed methods 
(Base, MC, Phase, PhaseMC) by comparing these tra¬ 
jectories to the true phantom motion trajectory. Each 
CBCT projection dataset (656 images) constituted 656 
measurements of the marker position in phantom geom¬ 
etry. For each measurement, the error was given as the 
magnitude of the 3D vector difference between the cal¬ 
culated position and the true position. The root-mean- 
squared trajectory error was given for each trajectory as 
the square root of the mean of the squared errors for all 
measurements in that trajectory. Additionally, the fre¬ 
quency of large errors was computed by calculating the 
fraction of measurements with errors exceeding 3 mm. 

The accuracy of dosimetric calculations was computed 
using our phantom-simulated clinical scenario. Since the 
dosimetric results depend on the shape of the target, 32 
abdominal/lung tumor volumes were used to evaluate the 
accuracy of the proposed methods (mean volume, 25 ± 
23 cm^; range 2.0 - 80 cm^). The reference dose was cal¬ 
culated by performing this analysis using the exact tra¬ 
jectory, and the calculated dose used the reconstructed 


trajectories. The mean dose error was calculated for each 
volume as the mean dose difference between calculated 
dose and reference dose. The mean absolute dose er¬ 
ror was calculated as the mean of the magnitude of these 
dose differences across all voxel of the volume. Mean dose 
error gives a sense of the dosimetric trends of each tra¬ 
jectory reconstruction method, and can reveal whether 
a given method underestimates or overestimates dose. 
Mean absolute dose error denotes the overall accuracy of 
a method in determining accurate dose in each voxel (by 
preventing positive and negative errors from cancelling). 
In addition to these metrics, dose-volume histograms of 
each volume were computed for each method, and com¬ 
pared to the reference. 

III. RESULTS 

A. Accuracy of trajectory reconstruction 

The accuracy of the trajectory calculations are shown 
in Fig. 3; results were calculated without respiratory 
phase or Monte Carlo corrections (Base), with phase cor¬ 
rections only (Phase), with Monte Carlo corrections only 
(MC), and with both phase and Monte Carlo corrections 
(PhaseMC). Trajectory amplitude ranged from 0 to 2.5 
cm. 

With all methods, the mean position of the marker was 
determined with accuracy better than 0.1 mm. In Fig. 
3a, it can be seen that the Phase and PhaseMC methods 
reduced the root-mean-squared (RMS) trajectory errors 
by 50% compared to the Base and MC methods, with 
an average RMS error of 3.8 ± 1.1% of the trajectory 
amplitude. In Fig 3b, the Phase and PhaseMC meth¬ 
ods decreased the fraction of large errors (defined here 
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FIG. 4. Mean and mean absolute dose error in the reconstructed trajectories. Errors were calculated using 32 clinically-drawn 
abdominal/lung CTVs. Motion amplitude ranged from 0 to 2.5 cm, and results were calculated with a range of isotropic margins 
and for trajectories with added noise. Dosimetric calculations using phase-resolved methods were more accurate, with mean 
absolute error less than 0.5% with zero margin, and with error less than 1% in the high-noise calculation. As noise increases, 
expectation-based methods generate a trajectory which over-estimates the dose to the CTV. Monte Carlo-based trajectories 
under-estimate dose by an average of 0.2% of prescription dose, even as trajectory noise increases. In A) and B), the noise is 0 
mm and the margin is 4 mm. In C) and D), the amplitude is 2.5 mm and the noise is 0 mm. In E) and F), the amplitude is 
1.2 mm and the margin is 4 mm. 


as trajectory error > 3 mm) by an average of 80%. MC 
methods increased both the RMS trajectory error and 
percentage of large errors slightly. Fig. 3c shows his¬ 
tograms of trajectory errors in all projections across all 
acquired datasets (amplitudes ranging from 1 mm to 25 
mm in 1 mm increments). Trajectory error was less than 
1mm in 92% of projections (Phase and PhaseMC) and 
82% of projections (Base). 

B. Dosimetric Accuracy 

Fig. 4 shows the distribution of mean dose errors 
and mean absolute dose errors using trajectories recon¬ 
structed with the Base, MC, Phase, and PhaseMC meth¬ 
ods. Results are shown across variations in trajectory 
amplitude, PTV margin size, and amount of added tra¬ 
jectory noise. 

Dose errors increased sharply as the motion amplitude 
increased, although in the Phase and PhaseMC methods, 
absolute error never exceeded 0.5% of the prescription 
dose. In general, dose errors decreased as the margin size 
increased, mainly due to an overall reduction in the con¬ 
tribution of motion seen by having more margin around 
the CTV. However, it can also be seen that mean dose 
error was lower when using Phase and PhaseMC trajec¬ 


tory estimation methods. This is likely due to the lower 
RMS trajectory error demonstrated by these methods. 
As gaussian trajectory noise increased, the expectation- 
based methods Base and Phase tended to overestimate 
dose, while Monte Carlo methods MC and PhaseMC 
tended to underestimate it. While Monte Carlo methods 
prevented the overestimation of dose, they also slightly 
increased the overall absolute error (<0.1%). 

Fig. 5 shows the accuracy of each method in de- 
terming dosimetric parameters D90 and D95 for vary¬ 
ing amplitudes, PTV margin sizes, and trajectory noise. 
One can see that the error of Phase methods was lower 
than the error in methods which ignore respiratory phase. 
Expectation-based methods again overestimated the dose 
(relative to the true value) and Monte-Carlo methods un¬ 
derestimated it, but the Base and Phase methods had a 
lower absolute error. The greatest accuracy was seen in 
the Phase method, which was able to determine D90 and 
D95 of the CTV to within 0.5% relative to prescription. 


IV. DISCUSSION 

In this work, we have developed and tested methods 
to evaluate dose to a moving tumor using data from a 
standard clinical cone-beam CT. In particular, we have 
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FIG. 5. Accuracy of calculating clinical dose metrics D90 and D95 from reconstructed trajectories. Errors are calculated in 
32 clinically-drawn CTVs. Motion amplitude ranged from 0 to 2.5 cm, and results were calculated with a range of isotropic 
margins and for trajectories with added noise. Phase-resolved methods increased the accuracy of these dose metric calculations. 
As trajectory noise increased, the Base and Phase methods tended to overestimate D90 and D95, while MC methods prevented 
overestimation of dose. In A) and B), the noise is 0 mm and the margin is 4 mm. In C) and D), the amplitude is 2.5 mm and 
the noise is 0 mm. In E) and F), the amplitude is 1.2 mm and the margin is 4 mm. 


tailored these methods to the problem of calculating dose 
to highly mobile tumors implanted with hducial markers, 
such as lesions in the lung, pancreas, or liver’^'*'’^®. The 
trajectory and subsequent dose received by the target 
volume was determined with minimal error. Errors in 
mean position did not exceed 0.1 mm, and the RMS tra¬ 
jectory error averaged 3.8% of the total range of motion. 
Using the Phase method, the mean dose and mean ab¬ 
solute dose were calculated with over 99% accuracy, and 
D90/D95 were determined to within 0.5% of the true 
value. 

By incorporating information related to the respira¬ 
tory cycle, the accuracy of tumor trajectory reconstruc¬ 
tion and subsequent dosimetric analysis was improved. 
The 3D Gaussian PDF utilized in this study stems from 
the work of Papiez and Langer^'*, where it was assumed 
that any systematic position errors are corrected with 
daily image guidance, and that only random fluctuations 
remain. Thus the central limit theorem dictates that the 
distribution of the remaining random errors is Gaussian. 
However, this assumption is broken for tumors which un¬ 
dergo regular respiratory motion, and the distribution 
of the trajectory about the mean position is no longer 
well definied as Gaussian. By creating a unique Gaus¬ 
sian distribution for each portion of the respiratory cy¬ 
cle, this effect is minimized, and the accuracy of trajec¬ 
tory reconstruction is improved. This method does have 


limitations; it cannot be used for tumors which undergo 
non-cyclical motion, such as prostate tumors. It is also 
difficult to determine the respiratory phase in real-time. 

Given the prevalence of GBCT, this method may pro¬ 
vide a useful technique to clinically estimate the dose- 
of-the-day received by a mobile tumor, since over half of 
radiation oncology clinics perform some cone-beam GT 
imaging for localization^®. Our analysis is directed to¬ 
wards tumors implanted with fiducial markers, but also 
holds in cases without fiducials as long as the tumor lo¬ 
cation can be discerned in each projection (as has been 
demonstrated in lung tumors®^’®®). In theory, our method 
allows for very accurate trajectory/dose reconstruction, 
since data are obtained using a much more comprehensive 
motion profile than other techniques. A GBGT projec¬ 
tion dataset constitutes upwards of 600 samples of tu¬ 
mor position, whereas 4D GT typically uses only 10 res¬ 
piratory phases, and has been shown to not accurately 
represent daily intrafraction motion of the abdomen®®’®'*. 
As radiotherapy becomes more conformal, and hypofrac- 
tionation becomes more commonplace, it becomes crucial 
to understand and account for the effects of tumor mo¬ 
tion on dosimetry. This method could be used to design 
patient-specihc, optimized margins for treatment, or to 
verify the adequacy of the margins used. It could also be 
used to decide between motion management strategies 
on a patient-specihc basis (such as breath-hold, gating. 
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or free-breathing). It should be said that while we have 
validated the accuracy of our method in a simple clinical 
scenario, true dose-of-the-day reconstruction would re¬ 
quire a more sophisticated dosimetric approach. By con¬ 
structing a simplified clinical scenario which maximizes 
the impact of tumor motion, our results indicate that 
tumor trajectories can be derived from CBCT projec¬ 
tion datasets with sufficient accuracy to determine true 
dose. An additional benefit of this method is that the 
processing time is short. Template-matching to deter¬ 
mine seed locations can be performed in real-time, and it 
takes roughly 2 seconds to calculate the trajectory once 
the projected fiducial locations are known. However, to 
build the 3D PDF, one must know the entire set of pro¬ 
jected locations; in other words, the PDF cannot be built 
until after CBCT acquisition is completed. Our method 
could be applied to generate an accurate estimate of tu¬ 
mor dosimetry shortly after the completion of the CBCT 
acquisition. 

The conditions of this study were chosen to test the 
boundaries of dose reconstruction in highly mobile tu¬ 
mors. Ranging up to 2.5 cm, the amplitude of motion 
in the SI direction corresponded to the displacement of 
very highly mobile lower lobe lung tumors, and exceeds 
the range of liver and other abdominal tumors’^^’^'^. Hys¬ 
teresis was also exaggerated, and motion in the AP and 
LR directions (up to 1.2 and 0.6 cm) exceeded that seen 
in clinical cases“’'*'*. In addition to the large range of 
motion, the effects of that motion on the dose were ex¬ 
aggerated by assuming perfect dose fall-off outside the 
PTV. Our methods were also tested over a wide range 
of tumor volumes, and with up to 1 mm of Gaussian 
noise added to each point in the trajectory. Even un¬ 
der these circumstances, dose metrics of the tumor were 
determined accurately. One may note that the RMS tra¬ 
jectory error reported in this work (1-2 mm) exceeds the 
mean error in the similar work of Poulsen et aP®. This 
can be explained by the large amplitude of motion in the 
current study, and agrees well with the maximum RMS 
error (> 2 mm) reported for lung tumors in that work. 
Additionally, these results give insight as to the scale of 
motion which causes the greatest dosimetric inaccuracy 
for each method. In Fig. 4C, errors increase for the Base 
and MC methods as the magin size approaches the range 
of 5-8 mm, while the Phase and PhaseMC see the great¬ 
est error for 4-5 mm margins. Since our dosimetric model 
is sensitive to motion equal to or greater than the PTV 
margin, and the frequency of motion errors decreases ex¬ 
ponentially for larger errors (Fig 3C), these results sug¬ 
gest that motion on the order of 7 mm has the grestest 
effect for the Base and MC methods, and on the order of 
4 mm for the Phase and PhaseMC methods. 

It is clear from these results that incorporation of respi¬ 
ratory phase information is beneficial when reconstruct¬ 
ing the tumor trajectory of highly mobile tumors. How¬ 
ever, the role of MC sampling is mixed. In the Base 
scenario, the position of the tumor is determined as the 
expectation value of the Gaussian PDF along the detec¬ 


tion ray. Yet one expects this determination to exhibit 
a Gaussian error profile, as the tumor position under¬ 
goes random fluctuations in addition to the cyclic respi¬ 
ratory motion. MC sampling addresses this by sampling 
a greater portion of the tails of this distribution (rela¬ 
tive to using the expectation value), thereby generating 
a distribution of positions with a more realistic propor¬ 
tion of large deviations. One may expect the MC method 
to increase trajectory error; however, under certain con¬ 
ditions one would expect a more complete sampling of 
the error distribution to result in a more accurate calcu¬ 
lation of derived metrics (such as target dosimetry). As 
the trajectory noise increases, the Base and Phase meth¬ 
ods to overestimate the true tumor dose, as the expecta¬ 
tion value of the distribution likely lies within the PTV. 
MC sampling slightly increases error in an absolute sense 
(Fig. 4), but is a conservative estimate of D90/D95 target 
dose metrics in noisy trajectories (Fig. 5). MC methods 
may be useful for calculating dose to the prescription vol¬ 
ume, and could constitute a worst-case scenario estimate 
of tumor dose. However, for organs-at-risk and other 
volumes, one would likely use expectation-based position 
estimation (i.e. Base or Phase), which yields a more ac¬ 
curate result. 


V. CONCLUSIONS 

We have developed methods for accurately determining 
the trajectory of highly-mobile abdominal tumors using 
cone-beam CT projections. These methods have been 
applied to reconstruct the 3D trajectory of a respiratory 
motion phantom, and to calculate the dose received un¬ 
der that trajectory. Results were calculated without res¬ 
piratory phase or Monte Carlo corrections (Base), with 
phase corrections only (Phase), with Monte Carlo correc¬ 
tions only (MC), and with both phase and Monte Carlo 
corrections (PhaseMC). In all methods, errors in mean 
position did not exceed 0.1 mm, and the RMS trajectory 
error was only 3.8% of the total range of motion. Using 
phase-resolved methods Phase and PhaseMC, the mean 
dose and mean absolute dose were calculated with over 
99% accuracy, and D90/D95 were determined to within 
0.5% relative to the prescription dose. These methods 
can be applied to estimate the dose-of-the-day for mobile 
tumors, and can aid in the selection of motion manage¬ 
ment strategies. 
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